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A SET OF FORTRAN IV ROUTINES USED TO CALCULATE THE 
MASS FLOW RATE OF NATURAL GAS THROUGH NOZZLES 
by Robert C. Johnson 
Lewis Research Center 

SUMMARY 

A set of FORTRAN IV subroutines is presented that calculates the isentropic mass 
flow rate of natural gas through nozzles and also thermodynamic functions such as com- 
pressibility factor, entropy, enthalpy, and specific heat. The pressure range for these 
routines is 0. 1 to lOOxlO 5 newtons per square meter, and the temperature range is 
from 200 to 400 Kelvin. Three sets of independent variables are permitted. In addition 
to the plenum pressure and plenum temperature, the other independent variable may be 
either the nozzle -exit pressure, the nozzle -exit Mach number, or the nozzle -exit tem- 
perature. 


INTRODUCTION 

When nozzles are used for measuring the mass flow rate of natural gas, the con- 
ventional isentropic flow equations do not apply. These equations only apply to a perfect 
gas. A perfect gas is defined as one whose compressibility factor is invariant, with a 
value of 1, and whose specific heat is invariant. This ideal condition is closely approxi- 
mated by a gas such as air at room temperature and at pressures up to a few atmo- 
spheres. Natural gas, on the other hand, cannot be considered perfect even at atmo- 
spheric pressure because of the specific-heat variation with temperature. At higher 
pressures, the compressibility factor variation also becomes important. 

Since natural gas is being considered as a fuel for aircraft, as well as for other 
propulsion and power systems, an accurate method for making mass flow rate calcula- 
tions is necessary for gas metering. 

In reference 1, isentropic flow calculations were made for natural gas using the 
Benedict -Webb-Rubin state equation (refs. 2 to 4). The result of these calculations was 
a critical-flow factor. By using this factor, the isentropic mass flow rate of natural 



gas through critical -flow nozzles can be calculated. A critical -flow nozzle is one that 
operates with a throat Mach number of 1. The output of FORTRAN IV routines used to 
make the calculations in reference 1 includes more than the isentropic mass flow rate 
of natural gas. In fact, a set of thermodynamic point functions are evaluated at both the 
plenum and the nozzle exit. These functions include compressibility factor, specific 
heat, isentropic exponent, and specific-heat ratio. 

In all these calculations, the plenum independent variables are pressure and tem- 
perature. The nozzle-exit independent variable is either pressure, Mach number, or 
temperature. 

Because of the versatility of these routines as previously indicated, they are pre- 
sented in this report. The routines operate over a temperature range from 199 to 401 K 
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and a pressure range from 0. 1 to 101x10 newtons per square meter. For a successful 
calculation, the components of the natural gas mixture have to remain in the gaseous 
state . 


BASIC EQUATIONS 

The routines in this report use three basic relations. The first describes the 
pres sure -temperature -density behavior of natural gas. This relation is referred to as 
the state equation. The second describes the ideal-gas, specific-heat variation with 
temperature. The third describes the saturated -vapor -pressure variation with tempera- 
ture. This last relation is used to determine that the fluid is a gas. These three rela- 
tions are discussed in detail in the following sections. 


State Equation 


Since natural gas is a mixture of many gases, it is useful to have a state equation 
whose coefficients can be determined from the coefficients of the individual components. 
Such an equation is that developed by Benedict, Webb, and Rubin in reference 2. This 
equation is as follows: 

2 2 ~ a i^ 

c a nP (1 + a 1 P )e 

p 5 + -1 - 1 (i) 



where 


2 



( 2 ) 


R = 


8314.41 


m 


m 2 / (sec 2 )(K) 


and 


N 

m =E x i ra i 

i=l 


(3) 


(All symbols are defined in appendix A. ) 

In reference 3, a set of rules is presented that permits the coefficients of equa- 
tion (1) to be determined for a gas mixture when the coefficients for the component gases 
are known. These mixing rules are 


a l = 
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E x i (a i,i m f )1/2 

i=l 
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In table I, the appropriate functions of the Benedict -Webb-Rubin coefficients and the 
molecular weight that permit easy use of the mixing rules (eqs. (4) to (11)) are presented 
for hydrocarbons with up to four carbon atoms, and for the dilutant gases, nitrogen and 
carbon dioxide. Hydrocarbons having more than four carbon atoms are not included. 

For pipeline natural gas, the presence of appreciable amounts of the higher hydrocarbons 
would be unusual. 


TABLE I. - FUNCTIONS OF BENEDICT -WEBB-RUBIN COEFFICIENTS FOR NATURAL GAS COMPONENTS 


Function 

Methane 

Ethane 

Propane 

Butane 

2 -Methyl 
propane 

Nitrogen 

Carbon 

dioxide 




Molecular weight , 

m 




16.043 

30.07 

44.097 

58. 124 

58. 124 

28.013 

44.01 

(a im V/2 

0.0774618 

0. 108631 

0. 148328 

0. 184396 

0. 184396 

0.08660497 

0. 1264947 

(a ? m)*' 3 

. 3492534 

. 3974298 

.459968 

.4991506 

.5162001 

. 3577881 

. 3667953 

(a 3 ni) 1,2 

4.754745 

7. 116558 

9. 140405 

11.0863 

11. 16732 

3, 81227 

5.79486 

(a 4 m) 1 ' 2 

524.4702 

1479.446 

2488. 837 

3478. 505 

3218.478 

267,9035 

1273. 766 

(a 5 mV 73 

. 1500773 

.2232212 

. 2823162 

. 3419966 

. 3488057 

. 1256056 

. 1324808 

(a 6 m 2 ) 1/3 

. 8444029 

1.614287 

2.260465 

2.841437 

2. 869004 

, 5662865 

. 926749 

(a 7 m 2 ) 1/3 

31.41978 

73.64101 

116.2798 

156.8145 

151.624 

18. 83293 

53. 5166 

(a 8 m 3 ) r3 

. 04991572 

.0624375 

.08454082 

. 1032721 

. 1024136 

.06631022 

.09234484 


The coefficients for the hydrocarbons were calculated from the Benedict -Webb- 
Rubin coefficients tabulated in references 2 and 4. The nitrogen coefficients were cal- 
culated from the Benedict -Webb -Rubin coefficients tabulated in reference 5. Benedict - 
Webb-Rubin coefficients for carbon dioxide were calculated by a least -squares technique 
using the compressibility data in reference 6. From these, the coefficients for table I 
were calculated. The coefficients aj to a g are consistent with temperature in Kelvin 
and density in kilograms per cubic meter. Because of the use of these density units 
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i cither th 3.11 the molar density units used in references 2, 4, and 5, molecular weight is 
involved in the mixing equations (eqs. (4) to (11)), while molecular weight is not involved 
in the mixing equations given in reference 3. Computationally, equations (1) to (11) are 
equivalent to those given in references 2 to 4. 


Ideal -Gas Specific Heat 

The ideal-gas specific heat of a natural gas mixture is represented by 


"v, ideal 
R 


" ^0 + 




The ideal -gas specific heat of the natural gas components is given by 


(12a) 



(12b) 


(An "ideal gas” is defined as a gas that has an invariant compressibility factor whose 
value is 1. However, unlike a perfect gas, the specific heat varies with temperature. 
Real gases approach the ideal-gas condition as the pressure assumes lower and lower 
values, providing no dissociation takes place.) 

For a natural-gas mixture, the ideal-gas specific heat can be expressed in terms of 
the ideal-gas specific heat of the components as follows: 



From this, it is evident that 
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i=l 


X A, i 


(13) 


(14) 
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TABLE II. - COEFFICIENTS FOR IDEAL-GAS. SPECIFIC-HEAT EQUATION FOR COMPONENTS OF NATURAL GAS 


Coefficient 


Methane 


Ethane 


Propane 


Butane 

2 -Methyl 
propane 

Nitrogen 

Carbon 

dioxide 

.0068 

-3.06092 

2. 50115 

2. 50447 

.60962 

6.08128 

-9. 72058xl0" 3 

-. 508557 

. 235295 

-. 593889 

1.03606xl0' 2 

.48403 

. 87536--10" 3 

1. 34513* 10" 2 

-4.43726X10" 3 

-3. 73057x10" 


1.07774X10" 2 

6. 8256X10" 4 

-2.52264x10-" 


-1.31759X10" 3 

0 

6. 14015x10" 


0 

0 

-4. 11664x10"' 


0 

0 

0 


'2 

'*4 

' J 5 

'*6 

'*7 


2. 79983 
.4285 
-.27518 
2. 58217x10' 
2.41658x10" 
-2. 51637x10" 
-8.24658x10" 
1.15233x10" 


-9. 85338 
19.6577 
-10. 1866 
1. 82674 
.246368 
-. 120205 
1.08075x10" 
0 


-16. 7968 
29.0846 
-13. 8109 
2.21984 
. 365514 
-. 15326 
1. 29667-10" 
0 


The values of i are presented in table II for the seven gases that comprise the com- 
ponents of natural gas. These values were obtained by a least-squares fit of tabulated 
data. These fits are valid for a temperature range of 200 to 400 K. The methane data 
are from reference 7. These data are the result of an analysis of spectroscopic data 
and agree with the data of the American Petroleum Institute (ref. 8) to within 1 percent. 
The ethane and propane data below 273 K are from reference 9. The ethane and propane 
data above 273 K are from reference 8, as are the butane and 2 -methyl propane data. 
The nitrogen data are from reference 10 and the carbon dioxide data are from refer- 
ence 6. 


Saturated Vapor Pressure 

The equation used to express the vapor pressure as a function of temperature is 



(15) 


The coefficients for ethane, propane, and butane were determined by a least -squares fit 
of the data in reference 9. The coefficients for 2 -methyl propane were determined by a 
least-squares fit of the data in reference 8, and the coefficients for carbon dioxide were 
determined by a least-squares fit of the data in reference 6. The vapor pressure of ni- 
trogen and methane need not be considered since the temperatures considered in these 
routines are all above critical. Table III lists the coefficients involved in equation (15). 
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TABLE III. - COEFFICIENTS FOR EQUATION FOR SATURATED VAPOR 


PRESSURE AS FUNCTION OF TEMPERATURE 


Coefficient 

Ethane 

Propane 

Butane 

2 -Methyl 
propane 

Carbon 

dioxide 

b 0 

-8. 76886 

-13. 83014 

-19.89223 

-10. 14642 

-65. 13333 

b l 

18. 78746 

16.45255 

18.41968 

8. 17872 

48.09596 

b 2 

-5.205866 

-.765418 

-.787275 

2.679815 

30.296025 

b 3 

. 538879 

-1.080231 

-.980618 

-. 944109 

-34. 13448 

b 4 

0 

.0642219 

-.0129045 

-.275245 

10.442646 

b 5 

0 

.0667237 

.0766147 

. 128236 

-1.071251 

b 6 

0 

-.0097026 

-.0094861 

-.0124255 

0 


In addition to equation (15), the routines require direct representation of tempera- 
ture in terms of pressure; that is, 



(16) 


Table IV lists the coefficients for equation (16). 


TABLE IV. - COEFFICIENTS FOR EQUATION FOR TEMPERATURE AS FUNCTION 
OF SATURATED VAPOR PRESSURE 


Coefficient 

Ethane 

Propane 

Butane 

2 -Methyl 

Carbon 





propane 

dioxide 

c 0 

C 1 

c 2 

c 3 

c 4 

°5 

0. 107196 
. 134425 
9. 57502x10 ' 3 
-1.26548X10" 3 
2. 14075x10'® 
2. 33577x10'® 

0. 729923 
9. 58931X10' 2 
6. 82915X10' 3 
-4. 18242X10' 4 
-3. 12624X10' 5 
3. 74086x10'® 

1. 51144 
-4.29761X10" 3 
7.92127X10' 3 
6. 85975X10' 4 
-1. 14435x10 " 4 
5.80657X10' 6 

-5. 888696 
2.473601 
-.291051 
1.49661X10' 2 
-2.41504X10' 4 
0 

7. 571723 
-1.022272 
4. 65321X10' 2 
0 
0 
0 


CALCULATION PROCEDURE 

For these calculations, six functions of the compressibility factor are used: 


Z (P,T) = Z = -P— 
[ PRT 


( 17 ) 
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(18) 


Z n (P,T) = Z + t(— ) =— ($2 


,3T/ p RP\3T/ p 


z ttt (p,t) = z +p/— \ = JL./iE\ 
111 \ap/ T rt \ 0p/ T 


Zj V (p 


; T) = J (Z n 


1) 


dp 


z v (p 


> T) = 


< z n - z i> 


dp 


Zvt (P,T)=tM = C T, Ideal - C v 

Vi 1 dT l p R 


(19) 


( 20 ) 


( 21 ) 


( 22 ) 


For an ideal gas, Zj, Zjj, and Zj n equal 1, and Zjy, Zy, and Zyj equal zero. 

In addition, two functions of the ideal-gas specific heat are used in these calcula- 
tions. These are 


7 

r(T) = = + y\/jL\ k + Kg 

1 / R T U \100/ Z-J k \100/ 15 

J k=l 


(23) 


s n m = j 


IXlMgM dT = 100 
R 


7 

I 

k=0 


P k / T \k+l 


k + 1 \100 / 


+ K. 


H 


In terms of and the ideal-gas entropy and enthalpy are given by 


^M= 5t (T) -ln/jp 
R 1 \RTJ 


(24) 


(25) 


H 


ideal 

R 


£ n (T) + T 


(26) 
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TABLE V. - VALUE OF INTEGRATION CON- 


STANTS USED IN EQUATIONS (23) AND (24) 


Component 

Constant Kg 

Constant Kjj 

Methane 

-2.42592233 

-794.255051 

Ethane 

-16. 722706 

-224. 353146 

Propane 

-24.4685144 

43.254680 

Butane 

-6. 81234692 

-859. 768636 

2 -Methyl propane 

-7.67222838 

-656. 575168 

Nitrogen 

-1. 20430845 

-699. 709835 

Carbon dioxide 

-. 54815092 

-702. 986595 


The terms Kg and K R in equations (23) and (24) are constants of integration for 
the indefinite temperature integrals in these equations. The constant Kg is chosen 
such that the ideal-gas entropy equals zero at a temperature of 200 K and a pressure of 
1x10 newtons per square meter. The constant K R is chosen such that the ideal-gas 
enthalpy equals zero at a temperature of 200 K. Table V lists the values of Kg and K 
for the components of natural gas. The values of Kg and K R for natural gas mixtures 
are given by 


7 

Kg = In m + E x i (K s,i- In nij) (27) 

i=l 

7 

K H=E X l K H,i < 28 > 

i=l 

In terms of these compressibility -factor and ideal-gas specific-heat functions, the 
following thermodynamic quantities can be expressed in terms of density and tempera- 
ture: 


| = Sj - in P - Z, v (29) 

| = «n * T < Z I - z v> ( 3 °) 

.K 
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_ ^v, ideal 
R R 



y = 



i 

z iii 




c 

v, ideal 
R 



k = 


P 

P 





(31) 


(32) 


(33) 


With minor modifications, equations (29) to (33) are from reference 11. The gas is 
assumed to flow from a plenum, where the gas is at rest, to the exit of the nozzle, where 
the gas velocity can either be subsonic or supersonic. This flow is assumed to be isen- 
tropic and one dimensional. The independent variables in the plenum are pressure and 
temperature. The independent variable at the nozzle exit can be either pressure, Mach 
number, or temperature. The first case to be considered is where the nozzle -exit in- 
dependent variable is temperature. 

Since the thermodynamic functions are in terms of density and temperature, it is 
first necessary to solve the following equation for the plenum density: 


Z 0 ( P 0 ' T 0* RT 0 

Since p Q is involved implicitly in equation (34), an iterative solution for p Q is nec- 
essary. This solution is as follows: 

First estimate of Pq: 


0,1 


RT r 


(35) 


Succeeding estimates: 


p 0,n “ P 0,n-1 + 


r dP 

9p, 


(P 


0 “ P 0,n-1 


(36) 
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where 


(aj) r RT o * z m( p o,n-i ,T o) 

When Pq n ^ and Pq differ by less than 1 part per million, the computation is con- 
sidered complete. 

The condition that the flow is isentropic leads to the following equation: 


S 0- S e 


= 0 = |,(T 0 ) - Ijfiy - in -2 - z iv (p 0 ,t 0 ) + Z IV (P e ,T e ) 

R p e 


( 38 ) 


Since this equation is implicit in p , an iterative solution for P g is necessary. This 
solution is as follows: 

First estimate of P g : 


in P e>1 = lnP 0 -« I (T 0 ) + € I (T e ) 


( 39 ) 


Succeeding estimates: 


In P 


e, n 


In P 


e, n-1 


+ 




where 


and 


3 (In P e ) 



1 


Z Il( p e,n-l’ T e) 



S ( p e,n-P T e) " 


( 40 ) 


( 41 ) 


( 42 ) 
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When the S e and Sq differ by less than lx 10"' , the computation is considered com- 
plete. 

For physically valid solutions of Pq and P g , it is necessary that Zp Zjp and 
Zjn be positive; that is, the density must be positive, the pressure must increase with 
temperature at constant density, and the pressure must increase with density at con- 
stant temperature. In addition, the partial pressure of each natural gas component has 
to be less than a constant times the saturation pressure of that component. The value of 
this constant is initially 1, but it can be set to other values under program control. 

At this point, the thermodynamic states at both the plenum and the nozzle exit are 
defined. Again, because of equation (38), the plenum and nozzle -exit entropy are equal. 
In order to calculate the mass flow rate per unit area at the nozzle exit, the nozzle-exit 
velocity has to be determined. This can be calculated by means of the following equa- 
tion: 


V e =<2 


H(P 0 , T o) 


H (P e > Te* 


| 1/2 


(43) 


The enthalpies in this equation are calculated from equation (30). The mass flow rate 
per unit area is then given by 


G„ 


P x V 
e e 


(44) 


If it is desired to calculate the Mach number at the nozzle exit, the speed of sound at the 
nozzle exit has to be calculated. The equation for the speed of sound is 


a 


\ HP ( 


’ T e) 


x Z 


[( p e ,T e) 


x R x T. 


1/2 


(45) 


Then the Mach number becomes 



(46) 


Other quantities such as enthalpy, specific heat, and specific-heat ratio are calculated 
at both the plenum and nozzle exit by using equations (30), (31), and (33). 

If the nozzle -exit independent variable is either pressure or Mach number, a nozzle- 
exit temperature has to be estimated such that the calculated nozzle-exit pressure or 
Mach number agrees with the prescribed nozzle-exit pressure or Mach number. To do 
this, a series of temperature estimates is usually required. Each estimate is used in 
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equations (39) to (42). This ensures that the flow is isentropic. The procedure for es- 
timating the nozzle-exit temperature for these two cases is now described. 

If the nozzle -exit pressure is the independent variable, the first estimate of the 
nozzle -exit temperature has to satisfy the following conditions: 


T e,l< T 0 


T e, 1 ^ T sat 


(47) 

(48) 


Condition (48) holds for each and every component of the natural gas mixture. These 
conditions take precedence over the following equation for the first estimate of the 
nozzle -exit temperature: 



(49) 


Equation (49) is what would result if natural gas were perfect, with y equal to 4/3. 
The second temperature estimate is given by 




x (p e - P e , 1> 


perf 


(50) 


where 



(51) 


Equation (51) is based on a y of 4/3. The other estimates are given by 


f T - T 

T = T , + | e » n -l e,n-2 » x , ) 

e,n e,n-l I ) * (p e p e,n-l' 

p e,n-l " p e,n-2 


(52) 


When the calculated nozzle -exit pressure agrees with the prescribed nozzle -exit pressure 
to within 1 part per million, the nozzle-exit temperature is considered to be determined. 
If the nozzle-exit Mach number is the independent variable, the first estimate of the 
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nozzle -exit temperature has to satisfy conditions (47) and (48). These conditions take 
precedence over the following equation for the first nozzle -exit temperature estimate: 



(53) 


Equation (53) is what would result if natural gas were perfect with y equal to 4/3. The 
second estimate is given by 


where 


T e,2 = T e,l + 



x (M 


M e,l) 


i perf 



Equation (55) is based on a y of 4/3. The other estimates are given by 


T = T + 

e,n e,n-l 


'T _ T 

e,n-l e,n- 2 

M e,n-l- M e,n-2 ; 


x ( M e - M n _f) 


(54) 


(55) 


(56) 


When the calculated nozzle-exit Mach number agrees with the prescribed nozzle-exit 
Mach number to within 1 part in 10 thousand, the nozzle-exit temperature is considered 
to be determined. 


RESULTS AND DISCUSSION 

Using the equations and iterative procedures given in the previous section, a set of 
computer routines written in the FORTRAN IV version 13 language were developed to 
calculate the isentropic flow properties that would result from flow through a nozzle. A 
brief description of these routines follows. These routines are identified by their deck 
name. Routines whose names start with the letters "NG” depend on the nature of nat- 
ural gas. The other routines are general and can be used for any real gas whose com- 
pressibility factor is a function of density and temperature. 
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Deck NGCOMP 


This subroutine calculates composition-dependent constants required by the other 
routines. These constants are 

(1) Gas constant, R (eq. (2)) 

(2) Molecular weight, m (eq. (3)) 

(3) Benedict -Webb -Rubin coefficients a 1 to a g (eqs. (4) to (11)) 

(4) Ideal-gas, specific-heat coefficients, /3 Q to (eq. (14)) 

(5) Integration constant for entropy, Kg (eq. (27)) 

(6) Integration constant for enthalpy, K H (eq. (28)) 


Deck RGASC 

This is the subroutine that includes the iteration procedures necessary to calculate 
the isentropic mass flow rate of natural gas through a nozzle. These procedures are 
represented by equations (34) to (56). The procedures are general and would apply for 
most real gases whose compressibility factor is given as a function of density and tem- 
perature. Naturally, it is necessary to provide the other routines with the appropriate 
basic equations. In other words, as far as this particular subroutine is concerned, it 
could be used for other gases. In addition to the mass flow rate per unit area, the’output 
of this subroutine includes entropy, enthalpy, specific heat, and compressibility factor. 


Deck RDATA 


This is a block data subprogram that supplies constants that mainly have to do with 
the convergence limits of the various iteration procedures in RGASC. These constants 
do not depend on the nature of the gas. 


Deck NGZETA 

The compressibility-factor functions Zj to Zyj , as defined by equations (17) to 
(22), are calculated in this subroutine. 
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Deck NGTEMP 


The nondimensional ideal-gas specific heat (C y idea p/R and the related functions 
and £jp as given by equations (12a), (23), and (24), are calculated in this routine. 


Deck NGLOG 

This is a logical function that tests whether the pressure and temperature lie within 
the range of the state equation (eq. (1)) and within the range of the ideal-gas, specific- 
heat equation (eq. (12a)). In addition, this routine also tests whether all the natural 
gas components are in the gaseous state. 


Deck NGTLG 

This routine, if necessary, will change the temperature such that it is above the 
condensation temperature of each and every natural gas component. 


Deck NGDATA 

This is a block data subprogram that supplies constants for the other routines. 

These constants depend on the nature of the gas. 

This set of routines requires 2137 storage locations, not including the library rou- 
tines. The execution time for a typical case on an IBM 7094-7044 direct couple com- 
puter is of the order of 0. 04 second. 

In appendix B, the rules for using these routines are described; and, in addition, 
the input and output computer variables are defined. The card listing of these routines 
is given in appendix C. 

As an example of the results obtainable by these routines, two figures are presented. 
These figures apply to a typical natural gas mixture. The mole fraction of the com- 
ponents of this mixture are 0.9272 methane, 0.0361 ethane, 0.0055 propane, 0.001 bu- 
tane, 0. 0007 2 -methyl propane, 0. 0218 nitrogen, and 0.0077 carbon dioxide. Figure 1 
illustrates the mass flow rate deviation from perfect -gas behavior. In order to compress 
the scale of the ordinate, the ratio of real-gas to perfect-gas mass flow rate per unit 
area is multiplied by the square root of the plenum compressibility factor. This assures 
that the value of this expression approaches 1 as the pressure drop across the nozzle 
approaches zero. This occurrence is independent of plenum pressure and temperature. 
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The expression for the perfect-gas mass flow rate per unit area is as follows: 


r 




(57) 


Equation (57) assumes a y of 4/3. Even at a pressure level of zero, where the com- 
pressibility factor is 1, the mass flow ratio is not constant. The reason for this is that 
as the nozzle -exit pressure is reduced, the gas speeds up and its temperature drops. 

This drop in temperature causes the ideal -gas specific heat to vary. This is not perfect - 
gas behavior and equation (57) does not apply. At a plenum pressure of 100x10^ newtons 
per square meter, the gas is neither perfect nor ideal. The value of the plenum com- 
pressibility factor under these conditions is 0. 8366. So, while the ordinate of figure 1 
is 1.017 at a Mach number of 1, the actual mass flow rate ratio G /G is equal to 

1 1 12 ^ ^ y 

In figure 2, two quantities are plotted. One is the specific-heat ratio, and the other 
is the isentropic exponent as defined by equation (33). At very low pressures these 
quantities are equal to each other. At the higher pressures they are not equal, but differ 
by as much as 12 percent. This illustrates that, for this case, errors of as much as 
6 percent would be involved if the specific-heat ratio rather than the isentropic exponent 
were used to calculate the speed of sound. 

It should be pointed out that although these routines calculate thermodynamic func- 
tions that, for the most part, are more accurate than the functions that would be calcu- 
lated using the perfect-gas assumption, errors are still present. The principal sources 
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Figure 2. - Specific-heat ratio y and isentropic exponent 
k as function of plenum pressure for a natural gas mix- 
ture. 


of these errors are inaccuracies in the Benedict -Webb -Rubin state equation. Some es- 
timate of these errors can be made for the case of pure methane. Reference 12 pre- 
sents a state equation for pure methane that is probably accurate to within 0. 1 percent 
over the range of pressures and temperatures covered by the routines in this report. 
Using this state equation, the compressibility factor, the specific heat, and the critical 
mass flow rate were calculated. The details of these calculations are beyond the scope 
of this report. A comparison of the highly accurate values obtained by the state equation 
in reference 12 with the values obtained by the Benedict -Webb -Rubin state equation used 
in this report is shown in table VI. Except for the low -temperature region, the com- 
pressibility factors differ by less than 0. 5 percent. The specific heats differ by as 
much as 12 percent at a pressure of 50x10"* newtons per square meter and a temperature 
of 200 K. However, for most of the range, the deviation is under 1 percent. The mass- 
flow-rate deviation seems to be under 1 percent over the permissible range of the cal- 
culations. It should be pointed out that these comparisons do not apply exactly for nat- 
ural gas mixtures. However, since natural gas consists primarily of methane, these 
comparisons should still apply approximately. 
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TABLE VI. - COMPARISON OF THERMODYNAMIC FUNCTIONS FOR METHANE CALCULATED BY 
EQUATION (1) WITH THESE SAME FUNCTIONS CALCULATED BY THE STATE 
EQUATION IN REFERENCE 12 a 


Plenum station 

Plenum station 

Compressibility 

z o - Z 0 

7 * 

Z 0 

Ratio of specific 

c - c* 
C* 

Mass -flow -rate 

* 

G - G 

temperature, 

pressure, 

factor 

heat to gas 

ratio at Mach 1 

G* 

T o- 

% 



constant 

p 


G* 


K 

, 9 . 

Z 0 

Z 0 





G 


N/ ni 


o 

o 

c pV r 


G perf 

G perf 


200 

50X10 5 

100 

0. 560 
. 367 

0. 549 
. 362 

0.020 

.014 

B 

14.66 
10. 54 


Conditions at the nozzle 
exit are such that con- 
densation could occur. 

250 

1" : 1 

0. 836 

0. 835 

0.001 

5. 51 

5. 39 

0.022 


1.098 

-0.003 



.689 

.686 

.004 

7.97 

7.94 


1.262 

1. 273 

-.009 

300 


0.918 

0.919 

-0.001 

4.96 

4. 93 

0.006 


B5I 



■ ! .1 

. 854 

. 855 

-.001 

5.79 

5. 75 

.007 


■ 


350 

50x10 5 

0.957 

0.958 


4. 96 

1 S3 

-0.002 

1.017 

1.020 

-0.003 


100 

.928 

. 930 

1 

5. 39 

■ 

-.004 

1.048 

1.050 

-.002 

400 

50X10 5 


0. 980 


5. 14 

B 

-0.004 



0 


100 

.966 

. 970 

-.004 

5.41 


-.007 

Hpllll 

P|§|§|f 

-.002 


a The asterisk denotes the quantity obtained by using the state equation in ref. 12. 


CONCLUDING REMARKS 

A set of computer routines are presented that provide a means of calculating the one- 
dimensional, isentropic mass flow rate of natural gas through a nozzle, and also such 
thermodynamic functions as the compressibility factor, specific heat, entropy, and 
enthalpy. The results indicate that the inaccuracy of these calculations is less than 
1 percent over most of the permissible pressure and temperature range. 

The design of the routines is such that it is easy to modify them for other gases. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, January 15, 1971, 

720-03. 
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APPENDIX A 


SYMBOLS 


3 ^ -| j 3 / 


j j > • 




t>Q > bj , . . . , 



G 

H 

K H 

*8 

M 

m 

N 

P 

p sat 


R 

S 

T 


T 


sat 


V 


x 


Z 

z r z ip . . 


• ,Z VI 


a 


coefficients for the Benedict -Webb-Rubin state equation (eq. (1)) 

coefficients for saturation -pressure -against -temperature equation 
(eq. (15)) 

p o 

specific heat at constant pressure, m /(sec )(K) 

specific heat at constant volume, m 2 /(sec 2 )(K) 

coefficients for saturation-temperature -against -pressure equation 
(eq. (16)) 

mass flow rate per unit area, kg/ (m 2 ) (sec) 

2 2 
enthalpy, m /sec 

constant of integration involved in equation (24), K 
constant of integration involved in equation (23) 

Mach number 
molecular weight 

number of components in natural gas mixture 

p 

pressure, N/m 

minimum pressure at which condensation takes place at a given tern - 
perature, N/m 2 

2 2 

gas constant, m /(sec )(K) 
entropy, m /(sec )(K) 
temperature, K 

maximum temperature at which condensation takes place for a given 
pressure, K 

velocity, m/sec 

mole fraction of a natural gas component 
compressibility factor 

functions of the compressibility factor as defined by equations (17) to 

( 22 ) 

speed of sound, m/sec 
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^0 ’ ^ 1 ’ • • • ’ ^7 
7 

k’ l n 


P 

Subscripts: 

e 

i 

ideal 

j 

k 


n 

perf 

0 


coefficients for ideal-gas, specific-heat equation (eq. (12a)) 
specific-heat ratio 

functions of ideal -gas specific heat as defined by equations (23) and 
(24) 

O 

density, kg/m 

nozzle-exit station 

species 

ideal gas 

species 

running index 
th 

n estimate in iteration procedure 
perfect gas 
plenum station 
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APPENDIX B 


DESCRIPTION OF NATURAL GAS ROUTINES 

Since natural gas is a gas mixture whose composition is not fixed, the first sub- 
routine that is referenced is one that calculates a set of composition -dependent con- 
stants for use in the other routines. For a given composition, this only has to be called 
once in a given run. The following statement references this routine: 

CALL BDATA(X) 

The subroutine used to calculate the thermodynamic properties of natural gas is 
referenced by the following statement. 

CALL RGAS (KK, PA,TA, AM, PB, TB, FLOW, KODE) 

For a valid computation, the following conditions have to be satisfied: 

(1) 199 K <T < 401 K 

(2) 0. 1 N/m 2 < p < 101X10 5 N/m 2 

(3) The partial pressure of any component has to be less than a constant times the 
saturation pressure. Unless specified otherwise, this constant will have the value of 

1 . 

Some of the variables in these programs are entered or returned through labeled 
common. Therefore, the following common statements should be in the main program. 

COMMON/ LD ATA/ XKV, R, XMW 

COMMON/LIMIT/EDA, EDB, ETP, ETM 

COMMON/OUT PUT/OUT (9), CONV(4), ZA(6), ZB(6), KODl(5) 

The following symbols apply to these routines: 

X This is a seven-element array. The elements in this array are propor- 

tional to the mole fraction of the natural gas components. The order in 
which the gases appear are as follows: (1) methane, (2) ethane, (3) pro- 
pane, (4) butane, (5) 2 -methyl propane, (6) nitrogen, and (7) carbon dioxide. 

KK Controls entry to, and exit from, the subroutine. If KK = 0, just the plenum 

properties are calculated. If KK = 2, both the plenum and the nozzle -exit 
properties are calculated. If KK = 1, just the nozzle -exit properties are 
calculated. For a given set of plenum conditions, at least one reference 
has to be made for KK = 0 or KK = 2 before a reference can be made for 
KK = 1. 

p 

PA Plenum pressure, p 0 , N/m 
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TA 

AM 

PB 

TB 

FLOW 

KODE 


XKV 

R 

XMW 

EDA 


EDB 

ETP 


ETM 


OUT(l) 


Plenum temperature, Tq, K 

Nozzle -exit Mach number, M 

’ e 

O 

Nozzle-exit pressure, p g , N/m 

Nozzle -exit temperature, T , K 

Nozzle-exit mass flow rate per unit area, G , kg/(m 2 )(sec) 

Indicates the independent variables to the subroutine. If KODE = 1, these 
variables are PA, TA, and PB. If KODE = 2, these variables are PA, 
TA, and AM. If KODE = 3, these variables are PA, TA, and TB. 

Constant referred to in condition (3). Unless specified otherwise, the value 
of this constant is 1. 

Gas constant, R, m 2 /(sec 2 )(K) 

Molecular weight, m 

Maximum fractional difference permitted between calculated plenum pres- 
sure and prescribed plenum pressure. Unless otherwise specified, EDA 
equals IxlCT®. 

Maximum difference permitted between the nozzle -exit entropy and the plen- 
um entropy. Unless otherwise specified, EDB equals 1X10 -6 . 

This is pertinent when the nozzle-exit independent variable is pressure, and 
represents the maximum fractional difference permitted between the cal- 
culated nozzle-exit pressure and the prescribed nozzle-exit pressure. Un- 
less otherwise specified, ETP equals lxlO~*\ 

This is pertinent when the nozzle -exit independent variable is Mach number, 
and represents the maximum fractional difference permitted between the 
calculated nozzle -exit Mach number and the prescribed nozzle -exit Mach 
number. Unless otherwise specified, ETM equals lxlO -4 . 

Actual mass flow rate G g divided by the perfect -gas mass flow rate 
G e,perf‘ The P erfec t-g' a s mass flow rate is defined as follows: 


G f 
e,perf 


Pp 



for M„ £ 1 
e ' 


G e,perf = 6732 


P0 


Vrt, 


(57) 


for M = 1 


(Bl) 



OUT (2) 
OUT (3) 
OUT (4) 
OUT (5) 
OUT (6) 
OUT (7) 
OUT (8) 
OUT (9) 
CONV(l) 

CONY (2) 


CONV(3) 


CONV(4) 


ZA(1) 

ZA(2) 

ZA(3) 

ZA(4) 

ZA(5) 

ZA(6) 

ZB(1) 


Nozzle-exit specific heat. C /R 

P, e 

Nozzle -exit specific-heat ratio, 

Nozzle -exit isentropic exponent, k g 
Plenum enthalpy, Hq/R, K 
Plenum entropy, Sq/R 
Plenum specific heat, C n /R 
Plenum specific-heat ratio, 

Plenum isentropic exponent, kg 

Plenum pressure as calculated by the plenum density. This pressure should 
equal PA. 

For KODE = 1, this is the nozzle -exit pressure as calculated from the 
nozzle-exit density. This pressure should equal PB. For KODE = 2, this 
is the nozzle -exit Mach number as calculated from the nozzle -exit velocity 
and speed of sound. CONV(2) should equal AM. For KODE = 3, CONV(2) 
is equal to zero. 

Indicates the degree of convergence achieved in calculating the plenum den- 
sity 


CONV(3) = 


PO 


z o p o RT o 


Indicates the degree to which the nozzle -exit entropy equals the plenum 
entropy. 


CONV(4) = S e - S Q 


Z l( p 0’ V 

z n (p 0) t q ) 

Z III^ P 0’ T o) 

Zj V (Po’ T o) 

Zy( p 0» Tq) 
Zyi( p 0> Tq) 
Z l( p e- T e) 
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ZB(2) Z n (P e> T e ) 

ZE( 3 ) Z m (P e , T e ) 

ZB(4) Z iy (p e , T e ) 

ZB (5) Z y (P e , T e ) 

ZB(6) Z yi (P e , T e ) 

The following symbols represent integers that are used to indicate if the calculation 
is valid. If all the integers are zero, a valid calculation has been performed. If these 
are not zero, the conditions are as follows: 


KODl(l) 

KODl(2) 

KODl(3) 


KODl(4) 

KODl(5) 


Equals 1 if the plenum conditions are out of range in either pressure or tem- 
perature. A value of 1 terminates the calculation. 

Equals 1 if the nozzle -exit conditions are out of range in either pressure or 
temperature. A value of 1 terminates the calculation. 

If KODE = 1, this quantity equals 1 if the nozzle-exit pressure fails to con- 
verge to the prescribed nozzle-exit pressure. If KODE = 2, this quantity 
equals 1 if the nozzle-exit Mach number fails to converge to the prescribed 
nozzle -exit Mach number. 

Equals 1 if the iteration procedure for the calculation of the plenum density 
fails to converge. 

Equal's 1 if the iteration procedure for the calculation of the nozzle -exit 
density fails to converge. 
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APPENDIX C 


CARD LISTING OF NATURAL GAS COMPUTER PROGRAM 


UBFTC Nil COM P 


THt FLU LOWING SUBROUTINE CALCULATES THE COMPOSITION DEPENDENT 
CONSTANTS THAT ARE USED BY THE OTHER ROUTINES. THIS SUBROUTINE HAS 
TC BE CALLED BEFORE A CALL r(I RGAS CAN BE MADE. 

SUBROUTINE BOAT A (X) 

DIMENSION XI 7), MOU 7 ) ♦ XM0LI7), SI 7), HI 7), BE (8), 6WR(8,7), CP ( 8 
1 , 7 ) 

COMMON / L PAT A/ XKV.R.MW 

COMMON /GDA T A/ RC , D 2 , GAMA, GAMB , G AMC , G AMD , GAME 
COMMON /PDATA/ F I 7 ) / Z D AT A/ l { 8 ) / T DAT A / A I 8 ) , H I , S I 
REAL MOL, Mw, MW 2 
EQUIVALENCE I All), BE I I) ) 

DATA MOL, XMOL/ 1 b.043 , 3G. 07, A A .097, 2*88. 124,28.013, A A. 01 ,2. 775 2 7262 
1 ,3. AO 3 52 8, 7 863 9175 , 2 *4,062 6 A 748 , 3. 3 3266869, 3. 78441688/ 

DATA S, I </ -2. 4259 22 33, -16.7 22 706,-24.4685144,-6.01234692,-7.67222 83 
18, - 1.2 0430845, -C .6 481 SO 92, -7 9 4. 255051,-224. 353146, 43. 2 546 8, -859,76 
28636,-656.575168,-699.709835,-702. 986595/ 

DATA CP/2 . 7998255,4. 2 849 981-1,-2. 7 51 805E- 1 , 2 . 582 17 1 IE-2, 2 . 4165 792E 
1-2,-2.5163 73 7 E- 3, -8. 2465805E-4, 1 „ 1 523272 E-4 , -4 . 85 33835 , 19. 65 7673 , - 
210. 186582, 1.8267443, 0.2463631 ,-0. 1202G4 8, L.080 7487 E- 2,0. 0,-16.7968 
307, 29.084569,-13 .81088 3, 2.2190327,0. 3655141,-v. 1 53 2602 , 1 . 29666 BE- 2 

4 .0. 0,-1.0067996, 4. 609 fcl 9,-0. 2352952, 4. 8753 55 8E-3,4*0. 0,-3. 0609193, 
56.0 812842,-0.5938891 , 1.345 13 IE-2, l.C 77738E-2, -1.31 759 33E-3, 2*0. 0,2 
6.501 145,-9. 720 581 E-3 , 1 . 0 360 5 6E-2 , - 4 . 437258E- 3, 6. 82 5 59 6E-4, 3*0. u, 2. 

75044 68 4, -0.5 68 5 567, 0.4 8403()2»-3. 730571E-2.-2. 522643E-2, 6.1401476E- 

83, -4.1 16635 7E-4 ,0.0/ 

DATA BWR/ .0774618, . 34925 34 , 4. 75474 5. 524 , 4702 , . 1500 773, .8444029,3 1 . 
1419 78, .049915 72, . 10863 1 ,. 39742 98, 7 . 1 16558 , 14 79 .446, .2232212, 1 .6142 
28 7, 73. 64101 , .0624375 , . 14 832 8 , . 4^996 8 , 9 . 1 40405 , 24 88 . 837 2 323 16 2, 2. 
3260465, 116.279 (3, .08454082,. 134 396, .4991506, 11. 0863, 3478. 505, .34199 
466, X. 841437, 1-56.8145, • 10 327 2 1, . 184396, .5162001 , 11 . 16732, 3218.4 78, . 

53438057.2.369004. 1 51 . 6 24 . . 1024 1 36 . .08660497. .3577881.3.81227.267.9 
6035 , . 1256056, . 566 286 5 , 18 . 8 32 93, . 066 3 102 2 , . 1264 947, . 366 7953,5.7948 6 

7.1273. 766. . 1324808 9 26749 . 5 3. 5 1 66 . . 09234484/ 

XX-2.0 
DO 1 N= 1 » 7 

1 XX=XX+X(N) 

DO 2 N = 1 , 7 
BE I I ) = 0 . 0 

2 F(N)=XIM)/XX 
PEI ; ) =0.0 

SI = ' .0 
HI =: .0 
M W = v . 0 
CG 3 N = 1 , 7 

S I = S I + E I N ) * I S I N ) - X MO L ( N ) ) 

)• I = : ! I + F IN) *h ( N ) 

MW = ”W<-F IN ) *MGL I N ) 

PEI 1 )=8E I 1 ) +F( M) *bWR I 1 ,N ) 
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o o no 


3 


DO 3 1=3,8 

BE ( I ) = BE ( [ )+F(N)*riWR ( I, '4) 

SI=SI+ALOG(MW> 

CO 4 N= 1 , 7 
DC 4 M= 1 » 7 

4 BE ( 2 > = BE t 2 ) + F ( N) *F (M > *( tUR( 2 ,N) +.BWR ( 2, M))**3 
MW2=MW**2 

Z( 1 ) =-RE( 1 ) **2/MW2 
Z ( 2 ) =B E { 2 )/( 8.0*MW) 

Z(3)=-BE(3)**2/MW 

Z(4)=-BE(4)**2/MW 

Z ( 5 ) =B E ( 5 ) **3/MtrJ2 

Z(6)=-BE(6)**3/MW2 

Z(7t=(BE(6)*BE{8))**3/(MW»MW?**2) 

Z(8)=BE(7)**3/MW2 

DO 5 N=l,8 

A ( N I = 0 » 0 

DO 8 M= 1 , 7 

5 A(N)=AlN)+F(M)*CP(N t M) 

R = 8 3 1 4 . 41/Mv, 

RC=3.45105E-5«MW 

RETURN 

END 


$ 1 8FTC RGASC 


THE THERMODYNAMIC PROPERTIES OF A NON IDEAL GAS ARE CALCULATED IN 
THIS SUBROUTINE. 

SUBROUTINE RG AS ( KK , P A A , T A A, AMM , f»B8 , TBB , FLOW , KO ) 

COMMON /OUTPUT/ OUT { 9 ) , CON V( 4 ) , ZA { 6 ) , Z B { 6 ) , KOD1 ( 5 ) 

COMMON / L D AT A/ XKV.R.XMW 
COMMON /LIMIT/ EDA, EDO, FTP, ETM 

COM MON /GO AT A/ RC, D2.GAMA, GAMft , GAMC » GAMP , GAME , GAME 
DOUBLE PRECISION CP , C S , Cn, CH A, CSA ,CS8 , CAB, LRHUA , LRHOB , DZ A , D ZB 
LOGICAL LGFN 
P A = PAA 
T A = T A A 
K KK = KK 
KODF =KO 

GO TO ( 1,2,3 ) , KODE 

1 PB= P8B 
GO TO 4 

2 AM= AMM 
GO 10 4 

3 TB= TBB 

4 IF (KKK.EQ.l) GO TO 17 
DO b N = 1 , 3 

Z A ( \ ) = 1 . 0 

5 Z B ( N ) = 1 . 0 
DU b N = 4 , 6 
Z A ( N ) = 0 • 0 

6 Z B I N ) = 0 . 0 
DO 7 N= 1 , 5 



7 K 0 D L ( M ) = 0 

DO d N=l,9 
B OUT ( N) =C . 0 

CC 9 M = 1 ,4 

9 C0NV(N)=0.0 

IF (LGFNI PA , TA, KODlt 1 ) , ZA) ) RETURN 
C 

C THE ITERATION PROCESS FOR CALCULATING THE PLENUM DENSITY FOLLOWS. 

C 

A = P A / ( RttTA ) 

R H 0 A = A 
K N = 0 

10 DO 13 MM= 1.50 

CALL ZETA ( i , RHOA, TA , Z A ) 

IF (ZA(3> . LL . 0 . 0 > GO TO 14 

CGN V ( 3 ) = 1 . 0— ( P A/ RHOA ) / ( Z A( 1)*R*TA) 

IF ( ARS ( CONV (3) ) .LT.EDA) GO TO 16 
AAA=( ZA ( 1 ) -A/ RHOA ) /Z A ( 3) 

11 IF ( 1 . 0 - A A A ) 12,12,13 

12 AAA=AAA/2 .0 
GO TO 11 

13 RHOA=RHOA*( 1.0-AAA) 

14 IF (KN. EQ. 1 ) GO TO 15 
RH0A = D2 *A 

KN= 1 

GO TO 10 

15 KOD 1 ( 4 ) = 1 

16 CALL ZETA ( 3 , RHOA , TA , Z A ) 

IF (LGFNIPA.TA, KODlt 1 ) ,Z A) ) RETURN 
C 

C THE PLENUM THERMODYNAMIC FUNCTIONS ARE CALCULATED OY THE FOLLOWING 
C STATEMENTS. 

C 

CV=CP( T A ) - ZA ( 6.) 

GA=ZA( 3 ) +ZA ( 2 ) **2/CV 
OUT (8) =GA/ZA ( 3 ) 

CUT (9) =GA/ZA{ 1 ) 

OUT ( 7 ) =C V*OUT ( 8 ) 

CHA-CH (TA)+D3LE(TA*( ZAtl )-ZA(5) ) ) 

CUT (5) =CHA 
C SA=CS ( TA ) 

LRHOA=DLOG( DBLE ( RHOA ) ) 

DZ A=DBLE ( ZA { 4 ) ) 

CUT (6) =CSA-LRHOA— DZA 
CONVtl ) = Z A ( 1 )*TA*R*RHCA 
C 

IF (KKK.EG.O) RETURN 

17 GC TO ( 18, 19,20) ,KODE 

18 A M = 0 . 0 
C 

C THE INITIAL ESTIMATE OF THE NOZZLE EXIT TEMPERATURE WHEN THE NOZZLE 
C EXIT PRESSURE IS GIVEN IS MADE BY THE FOLLOWING STATEMENTS. 

C 

T8=TA* ( PB/PA) **GAMA 
GO 10 21 
C 

C THE INITIAL ESTIMATE OF THE NUZZLE EXIT TEMPERATURE WHEN THE NOZZLE 
C EXIT MACH NUMBER IS GIVEN IS MADE BY THE FOLLOWING STATEMENTS. 

C 
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L9 TRAT-1 . 0+GAMB*AM»*2 
PB=PA/TRAT»»GAMF 
TB= I A/ TRAT 
C 

GO TO 21 

20 PB=PA* ( T3/TA ) **GAMF 
GO TO 22 

21 CALL TLOGIC ( P B * T B ) 

22 IF {TB.LT.TA.AND.PB.LT. PA) GO TU 23 
K0D1 (2 ) =1 

Return 

23 TB 1 =TB 
K0D1 ( 3 ) =0 
NN= 1 

DO 24 N=l,4 

24 OUT ( N ) =0 . 0 
CONV ( 2 ) =0 . 0 
C0NV14) =0.0 
DO 25 N=l,3 

25 Z B ( N ) = 1 . 0 
FL0W=0 . 0 

DO 26 N=4 » 6 

26 ZB (N >=0.0 

27 KOD 1 ( 5 > =0 

IF (NM.FQ. 1) GO TO 28 
IF (LGFN(PB,TB,K0D1(2),ZB) ) RETURN 
C 

C THE ITERATION PROCESS FOR CALCULATING THE NOZZLF EXIT DENSITY 
C FOLLOWS. 

C 

28 C SB=CS ( TB ) 

CAB=CSB-CSA+LRHOA+DZA 
LRH08=LRH0A+CSB-CSA 
DO 29 M=l,50 
RHOB=DEXP ( LRHOB ) 

CALL ZETA ( 2 , RHOB , TB , Z 8 ) 

DZ8=DBLE( Z8(4 ) ) 

CONV ( 4 ) =C A8-DZ B— LR HOB 
IF (ABS(C0NV(4) ) .LT.EDB) GO TO 30 

29 LRHOB=L RHOB + CONV ( 4 ) / Z B ( 2 ) 

KOD 1 { 5 ) =1 

30 IF (RHOA-RHOB) 31,31,32 

31 KOD 1 ( 2 ) =1 
RETURN 

32 CALL ZETA ( 3 , RHOB , TB , ZB ) 

C 

C THE THERMODYNAMIC FUNCTIONS AT THE NOZZLE EXIT CONDITIONS ARE 
C CALCULATED BY THE FOLLOWING STATEMENTS. 

C 

VV=? ,3D3*( CHA-CHI TB) -DBLLI TB*(ZB(1)-ZB(5)))) 

CV=CP( TB ) -ZB ( 6 ) 

GA=ZB(3)+Z8I2)**2/CV 
CUT (4) =GA/ ZB ( 1 ) 

C 

GO TO ( 33,37,41) ,KODE 

33 AM=ASQRT(VV/(ZB(1)*0UT(4)*TB)) 

IF (NN.NE.l) B1 = C0NV ( 2 ) 

CQNV(2)=RH0B*ZB(1)*R*TB 

PERR=PB/CONV(2)-1.0 



IF (ABStPFRR ) .LT.ETP ) GO TO 42 
IF (MVI.GT.2C) GC TO 3 6 
|\.N=vN+l 
C 

C THE SUCCEEDING ESTIMATES OF THE NOZZLE EXIT TEMPERATURE ARE MADE 
C BY THE FOLLOWING STATEMENTS FOR THE CASE OF A GIVEN NOZZLE EXIT 
C PRESSURE. 

C 

IF (NN-2) 35,34,35 
34 TB=TB*( 1.0 + G4MA*PkRR ) 

IF (TB.GE.TA) TB=0.999*TA 

TB2=TB 

GO TO 27 

3 5 TB=TBM TB2-TB1 ) * ( PB-CCNV (2 ) ) / (CONV( 2 )-Bl ) 

T B l =TB 2 
T82 = TB 
GO TO 27 
C 

36 KOD I ( 3 ) = 1 
GO TO 42 

37 p B= ZB ( 1 ) «TB*R*RHOB 

IF (MM.NE.l) B l = CONV ( 2 ) 

CONV(2)=ASQRT(VV/(ZB( 1 )*TB*01JT(4) ) ) 

IF (ABS(1.0-CGNV(2)/AM).LT.ETM) GO TO 42 
IF (NN.GT.20) GO TO 40 
NN=NN+1 
C 

C THE SUCCEEDING ESTIMATES OF THE NOZZLE EXIT TEMPERATURE ARE MADE 
C BY THE FOLLOWING STATEMENTS FOR THE CASE OF A GIVEN NOZZLE EXIT 
C MACH NUMBER. 

C 

IF (N.N-2) 39,38,39 

38 TB = TB* ( 1.0-GAMD*TB*AM*(AM-C0NV( 2) ) /TA) 

IF (TB.GE.TA) TB=0.999*TA 

TB2 =TB 
GO TO 27 

39 T8= IB* ( T32-T61 ) * ( AM-CONV (2 ) ) / ( CUNV ( 2) -61 ) 

TB I =TB2 

TB? =TB 
GG TO 27 
C 

40 KOD 1 { 3 ) = 1 
GO TO 42 

41 AM=ASQRT ( VVZ ( ZB( 1 )*OUT (4 )»TB ) ) 

P6 = ZB ( I ) *R*RHGB*TB 

CON V ( 2 )=0.C 

42 IF (L3FN(PB»TB,KGD1( 2), ZB) ) RETURN 
IF (VV.GT.O.O) GO TO 43 

KODI (2 ) =1 
RETURN 
C 

C THE ISENTROPIC FLOW PROPERTIES ARE CALCULATED BY THE FOLLOWING 
C STATEMENTS. 

C 

43 FLOw=PB*SQRT (VV/R)/(ZB(1)*TB) 

OUT ( 3) =CA/ZB ( 3 ) 

OUT (2) =CV*OUT ! 3) 

T8F=(PB/PA)*«GAMA 

IF ((AM.EQ.I.O) .AND. (K0DE.EQ.2) ) GO TO 4^ 
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noon 


FLOW I = P A » S CRT ( G A M E » ( PB/P A ) **GAMO ( 1. Q-TBF ) / ( R* TA ) i 
GO IQ 45 

44 FLOW I =P A* SORT ( RC/TA) 

45 GUT ( 1 ) = F L 0 to / F L 0 « I 
C 

GO TO (44,47,48) , KODE 

46 

TQH=TB 

RETURN) 

47 PBB=PB 
TBB=TB 
RETURN) 

48 AMM=AM 
P88 = PB 
RETURN 
END 


ilBFTC RDATA 

BLOCK DATA 

COMMON /LUATA/ XKV,R,XMW 

COMMON /LIMIT/ FDA,E08,ETP,ETM 

DATA X K V, EDA , ED3 , E TP , E TM / 1 . 0 , 3 * 1 . OE-6 , I . C-E-4/ 

END 


SIBFTC NGZETA 

TFF FUNCTIONS OF THE COMPRESSIBILITY FACTOR ARE CALCULATED IN THE 
FOLLOWING SUBROUTINE. 

SUBROUTINE ZFTA ( KK, P A , T A, Z ) 

COMMON / Z D A T A / A, 81, 32,8 8, B4,B5,B6,R7 
DIMENSION Z ( 6 ) 

K = K '< 

P = P A 

t = r a 
T 3= r**3 
B7T3=B7/T3 
P2=P*P 
P4=!’2*P2 

P5= P4*p 

AP 2 = A* P 2 
PFXP = EXP ( AP2 ) 

IF (K.EO.2) GC TO I 
BT1=(BI+B2/T+B3/T3)*P 
BT2 = I B4 + B5/ I ) »P2 
BT5=B6*P5/T 
ZA= 1 .D+BT 1+BT2+BT5 
Z. R=B7T3*P?» ( 1.C-AP2) *PEXP 
Z ( 1 ) =Z A + ZB 
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0 0.00 


/ A = I .0+2.0*8Tl+3.C*8T2+6.0*6T5 
Z» = B7T3*P2*PEXP* ( 3.0* ( 1 . ' -AP2 ) -2 .0*AP 2**2) 
Z(3)=ZA+ZC 
IF (K.EQ.l) RETURN 
1 BTl=(Bl-2 .0*R3/T3) *P 

GT2 = P.A*P2 
Z A- 1 .QvBT 1 + JT2 

ZB = -2. 0*B7 T3*P2»PEXP* ( 1 . ^-AP2 ) 

Z ( 2 ) =Z A+ZB 
ZA=‘'T1+BT2/2.C 

ZB=-B7T3*( 2. C*(PEXP-1.0) /A-P2*PEXP) 

Z(A)=ZA+ZB 
IF (K.EQ.2) RETURN 

ZA=- (Q2/T+3.0*B3/T3> * P-0 .3* B 5* P 2/ T-C . 2*BA*P8 / T 

ZB=1.5*ZS 

Z( 5 )=ZA + ZB 

ZA=6.3*B3*P/T3 

ZB=-2.3*Z8 

Z ( 6 ) = Z A + Z B 

RETURN 

ENO 


$ I BF TC NGTENP 

THE ICFAL GAS SPFCIFIC HEAT AND FUNCTIONS OF THF IUEAL GAS SPECEFIC 
HEAT ARE CALCULATED BY THE FOLLOWING ROUTINE. 

DOUBLE PRECISION FUNCTION CP(T) 

COMMON / T 0 A T A / A (8), HI, SI 
DOUBLE PRECISION S.XN.CS.CH 
K = 1 

1 S=T/10G. 000 

GO TO ( 2 , A , 6 ) , K 

2 C P = A { 8 ) 

CO 3 N= 1 , 7 
roi = c-N 

3 CO = CP*S+A(.VN) 

RETURN 
ENTRY G S I T ) 

K = 2 

GO TO I 

A CP=A(8)/7.C 

DC 8 N= I , 6 
A N=8-N 
XN=7-N 

8 C 0 = C P * S + A ( N M ) / X N 

C P = C P * S + A (1) * AL OG ( S > + S I 

RETURN 

ENTRY CHIT) 

K = 3 

GO TO 1 

6 CP = M 8 ) /8.0 

00 7 N=I ,7 
> N=6- N 
X N = N N 
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n o o o o 


7 CP=CP*S+A INN) /XN 
cp=cp» r+H i 

RETURN 

F N D 


$ I BF TC N G L C G 

TFF FOLLOWING FUNCTION IS USED TC DETERMINE WHETHER OR NOT THE FLUID 
COMPONENTS ARE A GAS AND WHETHER OR NOT THE PRESSURE AND TEMPERATURE 
CF THE FLUID LIES WITHIN THE RANGE OF THE STATE EQUATION. 

LOGICAL FUNCTION LGF IN ( P , T , J , / ) 

COMMON /PDATA/ X ( 7 ) 

COMMON /ldata/ xkv ,r , xmw 
COMMON / P M T M / TM AX ( 5 ) » P M [ N ( 5 ) 

DIMENSION Z ( 6 ) » FIS) , A( 7, 5) 

DATA A/-8. 7688639, 18. 767^64, -5. 2058657, 0.53887d7,3«0.0?-13. 830142? 
116.452554,-0. 765 4184,-1. '802 311 , 6 . 422 1 948E- 2 , 6 .67237E -2 , -9 . 70262 1 E 
2-3, -19. 892234, 18 . 4 1968? , -0 . 7 87 2747 , - . 9806 1 78 ,- 1 . 2904525E-2 , 7 . 661 46 

367E-2, -9. 4861 348E-3, - 10. 14641 8 , 8. 1787204, 2. 679814 7, -.944109 4, -.2 75 

42447. . 128236,-1 . 2425 4 6 E- 2, -65. 133331 , 48. 095955 , 30. 296025,-34. 1 3448 
5 ,10.442 646,-1.071251 ,0.0/ 

S=T/1D0.0 
J = 1 

L G F N = . TRUE. 

IF ( P. GT. 1 01 . 0E5 . OR. P . LT .0 . 1 . OR . S . L T . 1 . 99.0R . S . GT . 4.01 .OR . Z I 1 ) . LE . 

1 0 . 0 . 0R . Z ( 2 ) . L E . 0 . 0 .0 R . Z ( 3 ) . L E . C . 0 ) RETURN 
F ( 1 ) =X ( 2 ) 

F I 2 ) =X I 3 ! 

F ( 3 ) =X I 4 ) 

F { 4 ) = X ( 5 ) 

F ( 5 ) =X I 7 ) 

DO 2 1=1,5 
PX = P*F I I ) 

IF (PX.LT.PMINI I >.OR.S.GT.TMAX( I ) ) GO TO 2 
PLOG = A( 7, I ) 

DO L N = 1 * 6 
v= /-N 

1 PLOG = PLOG*S + A(M, I ) 

IF I PX.GT.XkV*FXP( PLOG) ) RETURN 

2 CONTINUE 
J =0 

LGEN=. FALSE. 

RETURN 

END 
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$ [ B F TC NC-TLG 

THIS SI. ROUT I \T WILL, IF NECESSARY, CHANGE THE VALUE OF- TA SO THAT IT 
IS ABC VC THE CCNDENSAT ION TEMPERATURE OF EACH AND EVERY FLUID 
C C M P C N F L T . 

SUBROUTINE T LOG I C (PA»TA) 

COMMON / P D A T A / XI 7) 

COM YON / P w T M / T M A X I 5 ) , P M I N I 5 ) 

DIMENSION FIS), AI6.5), SSI 5) 

DATA A/. 107). 95 8, . 134 A? 49,9. 5750 16703,-1.26547 72E-3,2 . 1407483E-5 ,2 
1 . 335769P-6, . 7299231 , 9.5893073E-2 , 6 . 829 1 5 23 E- 3 » -4 . 1 824 175 b-4 , -3 . 1 26 
2 24 3 IE- 5, 3. 74F8577E-6 , 1 . 5 1 1 439 6, -4 . 2976 141E-3 , 7 . 92 1 27 1 4E- 3 . 6 » 8 597 54 
32E-4.-1. 14 4349 8E-4.5 . 606 :,7 12 F-6 , -5 . eH869 -57 , 2 .47360 13 . 2 9105 1 , 1. 49 
466 1S9E -2, ~2. 4 15C 3 76E-4,0. 0,7.57 17232, -1.0222722, 4. 6532148E-2.3&0.0 
5/ 

P = P A 

S=TA/100.0 
FI 1 )=X (2) 

F I 2 ) =X C 3 ) 

F I 3 ) = X I 4 ! 

F ( 4 ) =X I 5 ) 

F ( 5 ) =X I 7 ) 

DO 3 N= 1 , 5 
PX=P»F ( N ) 

IF (PX.LT.PMININ).OR. S.GT.TMAX(N) ) GO TO 2 
PLX=AL0G IPX) 

S S I N ) = A ( 6 , N ) 

DO 1 M= 1 » 5 
vx=6-M 

1 SSI N ) = S S ( N ) * PL X+ A ( MX * N ) 

GO TO 3 

2 SSI 'J ) = S 

3 CONTINUE 

TA=100.0*AMAX1( SS(l) , SSI 2) ,SS( 3) ,SS(4) ,SS{ 5) ,S) 

IF I TA.LT. 199.0) TA=199.C1 

RETURN 

END 


1 1 8 FTC NGDATA 


BLOCK DATA 

COMMON / P M T M / TMAX (5 ) , PM INI 5 ) 

COMMON/GD AT A/KC » 02 , G AM A , 3AMB , GAMC , GAMD, GAME , GAMF 

DATA TMAX, PM IN/ 3. 055 5 , 3.6996,2*4.01, 3.04,2 l 8 3 DC. , 20 IOC . , 1981. , 39 C 6 

1.. 157800./ 

DATA GAMA , GAM6 , GAMC, G AMD , GAM E , GrtM F » D 2/ .25, . 16o666b6 7, 1.5, .33335 3 33 

13.8.0. 4.0.5.6/ 

END 
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